function FEA = FeasibilityExploitationArchive(FEA, Offspring, N, Add)
% The Feasible Exploitation Archive of CMOEA-CD
% Add = 1 --- environmental selection of SPEA2
% Add = 2 --- environmental selection of NSGA-II
% Add = 3 --- environmental selection of modified NSGA-III

%------------------------------- Copyright --------------------------------
% Copyright (c) 2025 BIMK Group. You are free to use the PlatEMO for
% research purposes. All publications which use this platform or any code
% in the platform should acknowledge the use of "PlatEMO" and reference "Ye
% Tian, Ran Cheng, Xingyi Zhang, and Yaochu Jin, PlatEMO: A MATLAB platform
% for evolutionary multi-objective optimization [educational forum], IEEE
% Computational Intelligence Magazine, 2017, 12(4): 73-87".
%--------------------------------------------------------------------------

% This function is written by Zhe Liu

    %% Dominance relation calculation
    FEA = [FEA, Offspring];
    NonDominated = DominationCal(FEA, 1);
    FEA = FEA(NonDominated);
    PopObj = FEA.objs;
    [P, M] = size(PopObj);  
    if P <= N
        FEA = FEA([1:P, unidrnd(P,[1, N - P])]);
    else
        PopCon = FEA.cons;
        Cons   = sum(max(0,PopCon),2);   
        if sum(Cons <= 0) <= N
            [~, index] = sort(Cons);
            remained_solution = index(1: N);
            FEA = FEA(remained_solution);
        %% Enviornmental selection
        elseif Add == 1
            FEA    = FEA(Cons <= 0); 
            PopObj = FEA.objs;
            [P, ~] = size(PopObj); 
            zmin   = min(PopObj, [], 1) - 1e-6;
            zmax   = max(PopObj, [], 1);        
            PopObj = (PopObj - zmin) ./ (zmax - zmin);
            Del    = Truncation(PopObj, P - N);
            remained_solution = 1 : P;
            remained_solution(Del) = [];
            FEA = FEA(remained_solution);
        elseif Add == 2
            FEA      = FEA(Cons <= 0); 
            PopObj   = FEA.objs;
            zmin     = min(PopObj, [], 1) - 1e-6;
            zmax     = max(PopObj, [], 1);        
            PopObj   = (PopObj - zmin) ./ (zmax - zmin);
            CrowdDis = CrowdingDistance(PopObj);
            [~,Rank] = sort(CrowdDis,'descend');  
            FEA = FEA(Rank(1:N));       
        elseif Add == 3
            FEA      = FEA(Cons <= 0); 
            PopObj   = FEA.objs; 
            [W, N]   = UniformPoint(N, M); 
            zmin     = min(PopObj, [], 1) - 1e-6;
            zmax     = max(PopObj, [], 1); 
            PopObj   = (PopObj - zmin) ./ (zmax - zmin);
            Distance = sqrt(sum(PopObj.^2, 2));
            Angle_Pop_to_W = sin(acos(1 - pdist2(W, PopObj, 'cosine')));
            S   = FEA;
            FEA = [];
            for i = 1 : N
                Fitness   = Distance' .* Angle_Pop_to_W(i,:);
                [~,index] = min(Fitness);
                FEA = [FEA,S(index)];
            end
        end
    end   
end

function Del = Truncation(PopObj,K)
% Select part of the solutions by truncation

    %% Truncation
    Distance = pdist2(PopObj,PopObj);
    Distance(logical(eye(length(Distance)))) = inf;
    Del = false(1,size(PopObj,1));
    while sum(Del) < K
        Remain   = find(~Del);
        Temp     = sort(Distance(Remain,Remain),2);
        [~,Rank] = sortrows(Temp);
        Del(Remain(Rank(1))) = true;
    end
end